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The dynamical chaos in Lennard- Jones toy models of heteropolymers is studied by molecular dy- 
namics simulations. It is shown that two nearby trajectories quickly diverge from each other if the 
heteropolymer corresponds to a random sequence. For good folders, on the other hand, two nearby 
trajectories may initially move apart but eventually they come together. Thus good folders are 
intrinsically non-chaotic. A choice of a distance of the initial conformation from the native state af- 
fects the way in which a separation between the twin trajectories behaves in time. This observation 
allows one to determine the size of a folding funnel in good folders. We study the energy landscapes 
of the toy models by determining the power spectra and fractal characteristics of the dependence 
of the potential energy on time. For good folders, folding and unfolding trajectories have distinctly 
different correlated behaviors at low frequencies. 



PACS Nos. 71.28.+d, 71.27.+a 
I. INTRODUCTION 

The notion of chaos in physical systems has several 
meanings. For instance, in the context of spin glasses 
it refers to the phenomenon of instability of the ground 
state against weak perturbations in the exchange cou- 
plings In the context of dynamical systems, on the 
other hand, it refers to an acute sensitivity of trajecto- 
ries to the initial conditions. Both of these meanings 
have relevance for understanding proteins and random 
heteropolymers . 

In sequences of aminoacids, the static chaos can be 
probed by investigating how mutations affect stability 
of the ground state. Several studies (^-^] have demon- 
strated that the native state of a random heteropolymer 
is unstable against mutations whereas there is stability in 
designed sequences provided the lengths of the sequences, 
N, are sufficiently small. It should be noted that the mu- 
tations may have strong effect on folding kinetics even in 
the case of short proteins || . Investigation of the effects 
of an aminoacid substitution in a protein is an essence of 
experimental procedures aimed at determining the tran- 
sition state in folding ||. 

Here, we focus on the notion of dynamical chaos in 
such systems and ask how do two folding trajectories re- 
late to each other as the system evolves from two nearby 
conformations. Furthermore, can such information pro- 
vide clues about the nature of the energy landscape of 
the system? 

We consider continuum space models, as opposed to 
lattice models, since the latter have intrinsically dis- 
cretized dynamics of a rather arbitrary nature. Specif- 
ically, we consider three toy ./V=16 off-lattice models 
that have been extensively characterized before and 
whose native states are shown in Figure 1. The Hamil- 
tonians of these systems arc defined in terms of the 



Lennard-Jones potentials (see Sec. 2) and the time evolu- 
tion is determined by Newton's equations that are solved 
by using the methods of molecular dynamics The 
first two of these systems, denoted in (7||] as G and R' , 
are two-dimensional whereas the third, denoted by H, 
is a Go-like model jn| of a three-dimensional helix ||. 
The quality of folding of these systems has been deter- 
mined by using thermodynamic fl^| and kinetic | fl3f cri- 
teria. Among the three systems, G is found to be a good 
folder, R' is a bad folder, and H has itermediate folding 
properties. Thus G is an analog of a protein and R' cor- 
responds to a typical random heteropolymer and H is a 
borderline case. The task of this paper is to compare the 
chaos related properties across this range of foldability. 

As a measure of the distance between two conforma- 
tions a and b we take S a b where 
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and where f* 1 is the position vector of the i'th monomer 
in the conformation a. This distance involves relative 
distances between the monomers and is bounded due 
the finite spacial extent of any conformation. Thus this 
distance cannot diverge as time grows which makes the 
characterization of the dynamical chaos unconventional 
in this case. 

Notice that once the system folds from a random con- 
formation to the native state, it may either stay in the 
immediate vicinity of the native state or it may depart 
further away and then keep visiting the native state. This 
translates into two possible scenarios for the asymptotic 
behavior of S a b'. either the distance saturates asymptoti- 
cally at a finite value - the case of H and R 1 - or it tends 
to zero - the case of G. Thus good folding properties 
are reflected in the small asymptotic values of 5 a b which 
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may be used as an alternative criterion for good foldabil- 
ity. The asymptotic saturation of S a b in bad folders is 
achieved at a much shorter time scale than one needed to 
establish the asymptotic tendency in good folder. Thus 
the bad folders can be said to be more chaotic than the 
good folders. 

Studies of the distance between two trajectories pro- 
vide information about the energy landscape available to 
a sequence. A complementary information can be ob- 
tained by studying individual trajectories. This is shown 
in Section 4 where we consider fractal and spectral prop- 
erties of the potential energy curve, E p (t), as seen on a 
trajectory as a function of time, t. The power spectrum, 
obtained by Fourier transforming E p (t), is found to in- 
dicate a correlated noise pattern and it shows sensitivity 
to a sequence in a way which is consistent with H being 
intermediate. Furthermore, the corresponding low fre- 
quency power law exponent is found to be markedly dif- 
ferent for folding and unfolding trajectories in systems of 
good foldability. The fractal dimensionality of the E p (i) 
curve is determined according to a procedure developed 
in We find that, for the well folding system, 

the temperature dependence of this fractal dimensional- 
ity has a dip around a temperature which is optimal for 
folding. No such dip arises in poor folders. 



II. THE SEQUENCES 

We start our discussion by defining models that we 
study. The sequences denoted by G and R' are two- 
dimensional versions of the model introduced by Iori, 
Marinari and Parisi fl6|| . Their native states are shown 
in Figure 1. The Hamiltonian is given by 
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where i and j range from 1 to N=16. The distance be- 



tween the beads, dij, is defined as \ fi 



fj | , where r. \ de- 



notes the position of bead i. d^ is measured in units of a 
the typical value of which is a = bA. The harmonic term 
in the Hamiltonian, with the spring constant k, couples 
the beads that are adjacent along the chain. The remain- 
ing terms represent the Lennard-Jones potential. In |jl6| , 
Aij = Aq + \/pr)ij, where Aq is constant and rjij's are 
Gaussian variables with zero mean and unit variance; /3 
controls the strength of the quenched disorder. The case 
of rjij — and Aq — C would correspond to a homopoly- 
mer with the standard Lennard-Jones interaction used 
in simulations of liquids. In Eq. (2.1) e is the typical 
Lennard-Jones energy parameter. We adopt the units in 
which (7=1 and consider k to be equal to 25e. Smaller 
values of k may violate the self-avoidance of the chain . 
The coupling constants Ay for system R' are listed in ref. 
§. These are shifted Gaussian-distributed numbers with 



the strongest attracting couplings assigned to the native 
contacts. For system G, A^ is taken to be 1 or for the 
native and non-native contacts respectively. System i?' 
has been shown to be structurally overconstrained and 
hard to fold. 
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FIG. 1. Native conformations of there sequences studied in 
this work. 

The helical system H has a native state, shown in Fig- 
ure 1, that mimics typical a-helix secondary structures. 
In this case the distances between beeds are assumed to 
have the length do = 3.8A. As one proceeds along the 
helix axis from one bead to another, the bead's azimuthal 
angle is rotated by 100° and the azimuthal length is dis- 
placed by 1.5 A. The Hamiltonian used to describe the 
helix is a Go-like modification of eq.2 and it reads |9(| 
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The first term is a backbone potential which includes the 
harmonic and anharmonic interactions 
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We take do = 3.8A, k\ = e and /c2 = 100c. The interac- 
tion between residues which form native contacts in the 
target conformation is chosen to be of the Lennard-Jones 
form 
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We choose so that each contact in the native struc- 
ture is stabilized at the minimum of the potential, i. 
e. aij = 2~ 1 / 6 d(j, where dfj is the length of the cor- 
responding native contact. Residues that not form the 
native contacts interact via a repulsive soft core poten- 
tial V NON , where 
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Here a = 2~ 1 / 6 d cut , d cut = 5.52. 

The time evolution of the sequences is determined by 
the fifth order predictor-corrector scheme JlIJ . The inte- 
gration step is chosen to be 0.005r, where r = ma 2 /e is 
the characteristic time unit and m is the mass of a bead. 
In order to simulate systems in contact with a heat bath 
of temperature T, we augment the equations of motion 
by the Langevin uncorrelated noise terms as described in 

mr = -Tr + F c + r] , (2.7) 

where F c — —\J r E p and 

(7?(0)77(t)) = 2Tk B T5(i), (2.8) 

where fee is the Boltzmann constant. We take T equal 
to 2. In the following, the temperature will be measured 
in the reduced units of e/ks- 
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FIG. 2. The temperature dependence of the folding times 
for sequences G, R' and H. The results are based on 100 start- 
ing random conformations. The arrows indicate the folding 
temperatures. 

Figure 2 shows the T-dependence of folding times for 
the three sequences studied. The folding time is de- 
fined as the median first passage time and the folding 
is declared to be accomplished if the distance to the na- 
tive state, S ns , becomes smaller than the native basin 
size S c . 5 C is defined by the shape distortion method 
U and it is equal to (0.2±0.01)ct and (0.09±0.01)ct 
for sequence G and R', respectively. For sequence H, 
S c = (0.6 ± 0.05)1 = (0.12 ± 0.01)cr 

The folding is the fastest at the temperature T m i n 
below which glassy kinetics set in. We determine that 
T mi „=0.15 ±0.02, 0.4 ±0.02, and 0.3 ±0.02 for G, R' 
and H respectively. Similar estimates of T m j„ were ob- 
tained with a Monte Carlo "dynamics" j7j (at a larger 
CPU cost). 



Socci and Onuchic Jl3J have proposed that what de- 
termines good foldability is whether the folding temper- 
ature, Tf, is outside of the range of temperatures where 
kinetics become glassy. Tf is defined as a temperature at 
which the equilibrium probability of being in the native 
state is 1/2. Here, we rephrase this criterion in the fol- 
lowing way: a sequence is a good folder if Tf is greater 
than T m i n , otherwise foldability is bad. We determine 
Tf through a Monte Carlo process and get values indi- 
cated by the arrows in Figure 2. The relative values of 
T m i n and Tf indicate what was announced in Section 1: 
G is a good folder, H is intermediate and R' is a bad 
folder. Characterization based on the specific heat and 
structural susceptibility yields a similar conclusion . 
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FIG. 3. The dependence of the probability, P, of falling 
into the native state as a result of quenching on the distance 
to the native state. The arrows indicate the values of 8f. 
The results are obtained by the Monte Carlo method and are 
averaged over 100 starting configurations. 

Another relevant property of a sequence is its charac- 
teristic funnel size, Sf. It can be estimated by generating 
random conformation of the system and then quenching 
them (by setting T = 0) and determining whether the 
resulting quenched state is native or not || . This allows 
one to estimate the probability, P, of getting to the na- 
tive state from a conformation which is 5 away form the 
native state. The critical value of S, S c above which this 
probability becomes smaller than one may be identified 
as a characteristic Sf. Our results on P are shown in Fig- 
ure 3 from which we can deduce that Sf ~ 0.55c, 0.35er 
and 0.16cr for G, R' and H, respectively (for H we took 
a = 5 A) . Thus the good folder G has a much larger fun- 
nel, as measured by the distance S nsi than the bad folder 
R' . A direct comparison of Sf for G and R' to that for H 
is not meaningful because the Hamiltonians and dimen- 
sionalities are different. Thus even though Sf /a for H is 
the smallest among the three sequences, its foldability is 
intermediate. 
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III. DYNAMICAL CHAOS 

In order to study the dynamic chaos we monitor two 
trajectories which evolve in time from two conformations 
which are initially separated by a distance AR. We chose 
AR = O.OOler. Smaller values of AR yield qualitatively 
similar results. The forces due to the Langevin noise are 
identical for both trajectories. 

Figure 4 illustrates what happens with two initially 
nearby trajectories as a function of time for systems 
G, R', and H respectively. The trajectories are char- 



acterized by the potential energy, radius of gyration, 
R g , and the number of established native contacts, N c . 
Two monomers i and j are declared to form a native 
contact if the distance between them is in the interval 
[0.9dfj, l.ldfj}. Independent of which of these quantities 
is used, we see that the two trajectories quickly come to- 
gether for system G but continue to be clearly distinct 
for a much longer time in the case of system R\ System 
H involves larger scales on the y-axis and, on a closer 
inspection, behaves like R\ 
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FIG. 4. Time evolution of the potential energy E p , the gyration radius R g , the distance to the native state 5 n a and the 
number of native contacts N c for two typical trajectories which have the initial departure AT? = 0.001. The left-hand, center, 
and right-hand parts of the figure correspond to sequences G, R' , and H respectively and all data refer to T = T min . The first 
and second trajectory are denoted by solid and dotted line respectively. The left arrow indicates when the last native contact 
appears for the first time, and the right arrow indicates when folding takes place, i.e. 8 ns < S c . 
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FIG. 5. The time dependence of 8 for G, H and R' at 
T — T m i„ for two nearby trajectories which start at two 
nearby unfoloded conformations. The solid portions of the 
curves are averaged over 1000 to 1500 pairs of such trajecto- 
ries. The dotted portions are averaged over 20 to 200 pairs. 
The circles indicate values of the median folding times. 



Figure 5 shows the time dependence of S at T = T min 
averaged over many pairs of starting trajectories. The 
time scales extend to much beyond the folding times. For 
7?' and H the distance 5 remains nonzero at these large 
time scales but for G — 5 — > 0. Thus for good folders all 
two initially separated trajectories eventually come to- 
gether. This is not so for bad and intermediate folders. 
In general, the better the folder, the longer it takes to 
establish the asymptotic behavior in 8. We can rephrase 
it by stating that i?' is more chaotic than H because the 
asymptotic large separation between the trajectories is 
established sooner. 

Figure 5 refers to trajectories which start in unfolded 
conformations. It is interesting to consider what hap- 
pens when the starting conformations come from a closer 
neighborhood of the native state. We generate such start- 
ing "points" by evolving a trajectory from an unfolded 
conformation to various stages chracterized by predefined 
values of S ns . At each state, we spawn a nearby compan- 
ion trajectory which is displaced by AR from the current 
conformation of the leading trajectory. The stages can 
alternatively be characterized by the numbers of estab- 
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lished native contacts but the 8 ns is more convenient to 
use when one deals with long heteropolymers. 

The first panel in the left-hand part of Figure 6 shows 
the time dependence of 5 for G at T — 0.1, 0.15, and 0.2 
for times significantly shorter than those corresponding 
to Figure 5. It is interesting to point out that the short 
time behavior of 8 a b may yield information about the size 
of the folding funnel. This can be achieved by studying 
two initially nearby trajectories which start at various 
locations in the folding funnel, i.e. at various distances 
5 ns away from the native state, as illustrated in Figure 
6. Notice that the initial placement affects the character 
of the initial evolution of 8. For 8 ns > 6^ s = 0.55ct and 
T = T min = 0.15 (the first panel of Figure 6), the two 
trajectories diverge. Otherwise, they reduce their rela- 
tive distance. This can be interpreted as a situation in 
which the twin trajectories are placed within the folding 
funnel. Thus the threshold value of S ns under the opti- 
mal folding conditions should be a measure of size of the 
folding funnel. In fact, this threshold value agrees with 
the funnel size as determined from Figure 3. For temper- 



atures above T m j„ the kinetic conditions deteriorate, the 
funnel desintegrates, and the threshold behavior in chaos 
disappears: the system becomes more chaotic. Below 
T m in, on the other hand, the funnel also fades away but 
the system is close to the quenching conditions. As the 
system evolves it becomes to be driven by energy mini- 
mization so the trajectories start to come closer together 
right away even for substantial values of S ns . 

The center and right-hand parts of Figure 6 show re- 
sults of a similar analysis for systems i?' and H respec- 
tively. For bad folders, even if the trajectories start ap- 
proaching each other at low T, as it happens for H (the 
bottom right-hand panel of Figure 6), they do not meet 
asymptotically because its is unlikely that they will be 
simultaneously in the same energy valleys. H displayes 
some borderline bahavior but there is no threshold be- 
havior in system i?' at any temperature: a viable folding 
funnel never forms. For good folders, however, studying 
twin trajectories allows one to establish the geometry of 
the folding funnel. 




FIG. 6. The time dependence of 5 for various types of starting configurations at temperatures which are indicated. The 
left-hand, center, and right-hand parts correspond to sequences G, R', and H respectively. For each T, the initial separation 
between two trajectories is AR — 0.001a. The solid line corresponds to the unfolded starting configuration. Sequence G: For 
T — T m i n = 0.15 the dotted lines correspond to 8 ns = 0.75, 0.6, 0.55, 0.5, 0.35 and 0.2a from top to bottom respectively. The 
curve with 5 n h 3 = 0.55a has a threshold character in the sense that for 8 ns > 8 n h 3 distance 8 decreases with time monotonically. 
The bar separates two different types of behavior of 8. For T = 0.2, the dotted lines correspond to 8 na — 0.7, 0.5, 0.3, and 
0.05a from top to bottom respectively. In the case of T — 0.1 the dotted lines are arranged in order S na = 1.6, 1.3, 1.1, 0.8 
and 0.5a. For T — 0.2 the distance between two trajectories grows for any S„ s . Sequence R': For T = T min = 0.4 the dotted 
lines correspond to S na = 0.8, 0.6, 0.4 and 0.05a from top to bottom respectively. For T = 0.2, the dotted lines are arranged in 
the order S na = 0.8, 0.6, 0.4, and 0.2a. For T=0.1 the order is S ns = 0.8, 0.7, 0.6 and 0.4a. Sequence H: For T = T min = 0.3 
the dotted lines correspond to 8 na — 2, 1.5, 1.2, 1, 0.8, 0.6, 0.4, 0.3 and 0.15a from top to bottom respectively. For T = 0.2 
the dotted lines correspond to 8 n3 = 1.2, 1., 0.8, 0.6, 0.5, 0.4, 0.2 and 0.1a. For T = 0.1 the dotted lines correspond to S na — 
1.5, 1.2, 1.0, 0.8, 0.6, 0.5, 0.4, and 0.2a. The results are averaged over 500 - 3000 starting configurations. 

In the case of R', for T = 0.2, the time dependence Namely, the apparent decrease of the distance is seen for 
of 6 is qualitatively the same as for T — T m j n = 0.4. <5^=0.4 a. It should be noted that the difference between 
The qualitative change is observed, however, at T = 0.1. good folder G and bad folder R' is clearly seen only at 
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T m i n . Away from T m i n they may behave qualitatively in 
the same way. 



IV. TIME EVOLUTION OF THE POTENTIAL 
ENERGY 

In order to understand the difference in the chaotic 
behavior of good and bad folders better we consider the 
time dependence of the potential energy, E p , on indi- 
vidual trajectories under the optimal folding conditions, 
i.e. at T m i n . We consider folding and unfolding trajecto- 
ries separately and demonstrate that the corresponding 
properties of E p are distinct. 
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FIG. 7. Spectrum of the potential energy for G, R' and H 
at T m in- The starting conformations are either unfolded or 
native as indicated at the top of the figure. The results are 
averaged over 400 trajectories. The values of / are shown 
next to the curves. 

The first question we ask is what are the properties of 
the power spectra of E p . Figure 7 shows the frequency 
dependence of \E p (uj)\ 2 , where 



(4.1) 



time scale of about 200r. The values of the power law ex- 
ponent are indicated in the Figure. For system G, there is 
also an intermediate time scale in the folding trajectories 
where s = 1.20 ± 0.20. It is possible that the existence 
of the intermediate frequency regime is a signature of a 
good foldability in general. 

The high frequency behavior corresponds to a large ex- 
ponent in the power law, of order 5-7 (the error bars 
here are of order 0.5). The low frequency i.e. long time 
behavior, however, is clearly distinct for G and almost 
the same for R\ For G, the low frequency behavior 
of E p corresponds to /=2.16 ± 0.20 and 0.68 ± 0.2 for 
the folding and unfolding trajectories respectively. For 
R\ on the other hand, one gets the 1// noise for both 
kinds of trajectories. This clearly indicates lack of any 
folding direction in R' . System H has an intermediate 
behavior again: for the low frequency folding trajecto- 
ries s = 1.88 ± 0.24 and for the unfolding trajectories 
s = 1.18 ± 0.20. Thus there is a difference between fold- 
ing and unfolding but the difference is not as strong as 
for G. 
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FIG. 8. Dependence of V(e) on e for G, R' and H. The re- 
sults are obtained at T = T m i n . The left and right hand pan- 
els are for folding and unfolding conformations respectively. 
The time steps are equal to the folding times at T = T m in- 
The values of 7 are shown next to the data points. 



The first observation is that if E p is viewed as noise then 
this noise is clearly correlated - there is no frequency 
regime which would correspond to white noise. In other 
words, if 



\E n 



(4.2) 



then the exponent s is non-zero. For each of the systems 
studied we observe two regimes with a power law depen- 
dence on lo: the low and high frequency regimes which 
are separated by luq ~ 0.032/r which corresponds to a 



Another way to characterize trajectories has been re- 
cently proposed by Lidar et al. |Q and it involves deter- 
mination of the fractal dimensionality, 7, which relates to 
the self-affinity properties of the E p (t) curve. This frac- 
tal dimensionality may be obtained by the e-variation 
method developed by Dubuc et al. p5[. Typically, for 
any function g one can introduce its e-variation ^(£,17) 
as follows 



V(e,g) = 



/ v(x, 
Jo 



e)dx 



(4.3) 
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where the e- variation is 



v(x,e) = sup g(x') 

R e (x) = {se[0,l];\x-s\<e} . (4.4) 



x'£R e (x) 



Figure 8 shows the e-dependence of e-variation V(e) 
for E p (t) of three typical trajectories at T m , n run in the 
folding or unfolding modes. In the folding mode, the tra- 
jectories are evolved until the folding is accomplished. In 
the unfolding mode, the trajectories are analysed for a 
duration of a typical median folding time at T m i n . The 
slopes in Figure 8 give the values of 7. If the starting 
conformation is native then 7 appears to be system inde- 
pendent and equal to about 1.90. This again indicates a 
correlated behavior since for Gaussian distributed num- 
bers 7=2. On the other hand, if the starting conforma- 
tion is unfolded, then 7 w 1.81, 1.94 and 1.93 for G, R' 
and H respectively. This suggests that the smaller the 7 
the better the foldability (and less chaos). 
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FIG. 9. Temperature dependence of 7 for G (open circles), 
R' (closed squares) and H (closed triangles). The starting 
conformations are unfolded. The time steps are equal to the 
folding times at T = T m i n . The number of starting configu- 
rations is 20 - 100. 



roughness. Therefore, 7 and s may depend on the sys- 
tem in a different way. It is interesting to ask why are 
the system-dependences of s and 7 in the folding mode 
stronger compared to the unfolding mode. The reason 
seems to be that the time dependence of the potential 
energy in the folding mode (shown in Figure 4) is sub- 
stantially stronger than in the unfolding mode for all of 
the three sequences. Thus, in the unfolding case the 
system-dependence becomes weaker and 7 even loses the 
dependence entirely. It would be interesting to know if 
this observation is still valid in real proteins. 

The temperature dependence of 7 for three sequences 
is shown in Figure 9 where the starting conformations are 
unfolded. The values of 7 presented in this figure fulfil 
the relation 7 = 2 — a, where a is the roughness expo- 
nent juj, if a is determined directly. Studies of models 
of real proteins H] indicate that 7 may depend on T 
weakly. This is also true for our model systems as shown 
in Figure 9. Note that system H again behaves in a way 
which is intermediate between i?' and G. Interestingly, in 
the protein-like sequence G we observe a dip in 7 around 
T m i n . A similar but wider dip was also observed for real 
proteins like myoglobin, BPTI and PPT juj. The pres- 
ence of the dip could be explained in the following way: 
around T m i n the system establishes a folding funnel and 
the motion becomes less rugged or less chaotic. Thus, 
fractal analysis around T m j n may provide the useful in- 
formation about the foldability. It should be noted that 
our results have been obtained for times equal to the fold- 
ing times at T min and the related conclusions are valid on 
these time scales. Longer or shorter runs may in princi- 
ple, change the estimates of 7 JlJ] . Within the error bars, 
the results shown in Figure 9, however, do not change if 
the the time scale is doubled. 

In conclusion, we have studied the dynamic chaos of 
several model sequences and demonstrated that good 
folders are essentially non-chaotic and bad folders are in- 
trinsically chaotic. The energy landscape of heteropoly- 
mers can be characterized by the spectral and fractal 
properties of the time evolution of the potential energy 
of the system. 

We thank T. X. Hoang for useful discussions and tech- 
nical help. This work was supported by Komitet Badan 
Naukowych (Poland; Grant number 2P03B-146 18). 



From Figures 7 and 8 one can see that s has a stronger 
system-dependence compared to 7. This is due to the fact 
that the spectral and roughness properties refer to differ- 
ent aspects of the behavior. This may be seen clearly 
in the case of the white noise where 7 = 2 but s = 
which means that the white noise is not correlated but 
its profile remains rough. In other words the spectral 
analysis provides information about the pattern corre- 
lation whereas the fractal dimensionality relates to the 
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